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Abstract 



In this series of papers we present a detailed study of the particle-particle 
collective excitations of the Hubbard model, and their contribution to the 
density and spin excitation spectrum. In the first paper, we shall investigate 
the singlet particle-particle pair with momentum (vr, tt), the rj particle, of the 
negative-C/ Hubbard model. We review three previously obtained theorems 
about the rj particle and develop a self-consistent linear response theory which 
takes into account its contribution to the density excitation spectrum in the 
superconducting state. We show that this self-consistent theory agrees with 
the exact theorems as well as the results of numerical Monte Carlo simulations. 
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I. INTRODUCTION 



The Hubbard model provides a basic framework within which the nature of strongly- 
correlated electron systems have been studied. Recently, it was found that the Hubbard 
model has 5*0(4) symmetryS, for both positive and negative on-site Coulomb energy U. 
That is, besides the usual SU{2) spin-rotation symmetry, it is also invariant under an SU{2) 
"pseudo-spin" rotation group which contains the charge U{1) symmetry as a subgroup. 
Based on this new symmetry, one of usna derived a number of theorems and argued that 
this symmetry implies the existence of a collective excitation, the 77-mode. This mode is 
characterized by a charge number of two, a spin quantum number zero, total momentum 
(tt, tt) and a sharp energy U — 2fi. Here /i the chemical potential. As the momentum shifts 
away from (tt, tt), the mode disperses and eventually merges into the continuum. Now, since 
this excitation has a charge quantum number of two, one can not couple to it in the normal 
state with the usual electric or magnetic-field probesQ used to study the density and spin 
excitations of a many body system. However, if the system were to become superconducting, 
spontaneously breaking the U{1) charge symmetry, then one should observe this rj mode in 
the charge density response. Thus one has the unusual situation, that above Tc there is a 
well defined excitation in the many-body system which is invisible to the usual probes and 
it is only below T^. that we see this excitation. 

Here we explore this phenomena for the negative-t/ Hubbard model. We begin with a 
review of the exact theorems regarding the ?7-mode. We then use a self-consistent linear 
response theory to approximately calculate correlations of the rj and the charge density p 
in the superconducting state. We show that this approach gives results in agreement with 
the exact theorems and hence provides a meaningful approximation. Finally, we carry out 
Monte Carlo simulations for the two-dimensional negative-f/ Hubbard model and compare 
these with the approximate solution. In the Monte Carlo case, we can calculate the 77-77 
correlation function above the Kosterlitz-Thouless temperature and clearly see the sharp rj- 
mode in the normal state in this correlation function. We can then calculate the p-p charge 
density response and see that the ?7-mode becomes visible in the charge density response 
when the temperature is lowered and the pairing correlations extend over the lattice. We 
discuss these results and their possible relevance to other problems in the conclusions. 

The formalism of self-consistent linear response theory takes into account the cou- 
pling of the particle-hole excitation and particle-particle excitation in the superconducting 
state. Such a formalism was first developed in the context of gaug^ invariance problem in 
super conduct ivity§^. It was later used by Bardasis and SchrieffeiB'l to discuss the possi- 
bility of excitonic states inside the BCS energy gap. While the formalism we use here is 
similar to these previous works, the physical origin of the rj mode is very different from the 
excitonic states inside the BCS energy gap. The excitonic states owe their existence to the 
BCS energy gap, they exist near a total momentum near zero, and disappear completely 
above T^. The 77 mode, on the other hand, owes its existence to the collapse of the particle 
particle continuum near (tt, tt) and exists even in the normal state, although they decouple 



In principle, a two-particle tunneling process could couple to this excitation but this would also 
require a superconductor. 
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from the density spectrum above Tc. After the prediction of the t] mode in the Hubbard 
modeli, Kostyrko and MicnaJl have apphed the conserving approximation in Nambu for- 
mahsm to study the collective modes of the extended Hubbard model at zero temperature, 
and numerically confirmed the existance of the t] mode in their formalism. In our present 
work we shall make the detailed comparison of 5*0(4) Ward identities and the correlation 
functions derived using the self-consistent linear response theory, invstigate the effect of 
finite temperature on the rj mode, and demonstrate explicitely the scaling of its intensity 
with the superconducting order parameter. 

Our present investigation is motivated in part by the recent neutron scattering exper- 
iments on the high superconductors, in which a collective resonance was found below 
the superconducttag trau^.tion te,nperatur<IHi^ This mode was interpreted by two of usij 
as a triplet particle particle collective mode, hereafter called the vr mode, of the positive 
U Hubbard model. Our motivation to study the details of the rj mode of the negative U 
Hubbard model is to use it as a theoretical laboratory to verify a mechanism in which a 
particle particle collective mode can couple to the particle hole spectrum below the super- 
conducting transition temperature, and to check the methodology used in the calculations. 
The 7] and vr modes both have charge two and exist near momentum (tt, tt), but they have 
different spin quantum numbers. The 7] mode has spin zero and energy of the order of \U\, 
while the vr mode has spin one and energy of the order of J, the spin exchange energy. The 
similarities and the differences of these two collective modes will be addressed in detail in 
our next paper. 



II. REVIEW OF THREE EXACT THEOREMS 

In this section we shall review three main theorems of the Hubbard model derived by 
one of us in referencesili. These results follow from Yang's wor k0 on the 7] pairing, and the 
5*0 (4) symmetry of the Hubbard mo delS. However, it is important to point out a physical 
difference between the following discussion and the idea of t] pairing. While the idea of rj 
pairing0'El refers to a ground state property of the model, the following discussion assumes a 
conventional, zero total momentum pairing in the ground state, and the t] pair is considered 
as an excitation of the system. 

Let us start with the Hubbard Hamiltonian defined by 

<ij> i i 

Within this model, one can define the following operators 



t _ \^ t t 

p 



?7 = {r]^ 

Vo = liNe~N) (2) 

where Q = (vr, vr), A^^e is the total number of electrons and A^ is the total number of lattice 
sites. It is easy to see that they form an SU (2) algebra, called the pseudospin algebra: 
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(3) 



These t] operators have a remarkable commutation property with the Hubbard Hamiltonian 

[H, r^t] = (f/ - 2fi)r]^ , [H, Vo] = , [H, r/^] = (4) 

Since the Hubbard Hamiltonian commutes with the Casimir operators t]q and 77^ = 
^{rj'^ri + riri'^) +770^70 of the pseudospin algebra, the Hubbard Hamiltonian posses an 50(4) = 
SU{2)spin X SU {2) pseudospin Symmetry group. These commutation rules and the symmetry 
property can be used to derive the following three important theorems, initially obtained 
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i. 



First let us consider the pure particle-particle correlation function involving only the r] 
operators: 



E 



(n|r/t|0) 



ijJ — EnQ + iO U + EnO — iO 



(5) 



where 77^ = Z]pCp+g|cLp| and the summation in the second line is over all the excited states 
|n). 

For the case of q = Q, the operator rjg reduces to the rj operator defined above. Through- 
out this paper we shall assume without loss of generality that system is less than half-filled, 
i.e. Ne < N. From the commutation relation with the Hubbard Hamiltonian, it is simple 
to see that r/"'' operator creates a single excited state when acting on the ground state,0 



\r]) = (1 - n)-^ r]^\0) 
where a normalization factor is introduced such that 
(r7|r7) = (l-n)-i(0|r7r/t|0) = (l 



(6) 



|0) 



(7) 



It is also easy to show that rj annihilates the vacuum: ri\0) = (we have already used 
this fact in (|^ to replace the product of rj and rj'' by their commutator). 

Therefore only one term survives in (^) with {n\ = {ri\. And for P{Q,uj) we have 



P{Q,u;) 



n] 



(1 — n) C(j — C(Jo + ^0 uo — ujq + i^i 



(8) 



^ In the consequent formulas we will be using the filling factor n = N^/N rather than A^e and 
will be omitting all the factors of for brevity. So, that summation over all momenta p in our 
notations gives 1 and not N. This may be viewed as analog of normilization by a unit volume in 
the continious case. When needed the explicit dependence on may be easily recovered. 
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where uq = U — 2fi is the energy of the r^-particle. Equation (H) is the content of our first 
theorem. 

Next let us consider a mixed correlation function involving the rj operator and the density 

operator Pq = Efca 



M{q,uj) = -iJ dte''^' < Tp^g{t)ril{0) > 

_^i {0\P~,\n){n\vl\0) (0|r^t|^)(^|p„^|o) -| 

For q = Q we have contribution from the |?7)-state only, and M{Q,uj) becomes 

MW..)^^"!''-'"'!"' ' . ^-2<^> ' (10) 



(1 — n) iu — luq + iO U{1 — n) iu — luq + iO 



where A = —UJ^k'^li'^-k] the superconducting order parameter. Equation (ITO) is our 



second theorem. 

Finally, let us consider the pure density-density correlation function defined by 



D{q,u) = -ij dte''^' < Tp^q{t)pq{0) > 



For D{Q,uj) we can pull out the part singular at luq 

D(Q,uj) = +part regular at ujq 

uj — coq + iO 



uj — ujq + iO 
< A >2 1 



+ part regular at cuq 



= '^TF?Ti T r~7T + P^"^^ regular at Uq (12) 

U''{1 — n) LJ — ujq + lO 

Equation ([l2D is our third theorem. 

From expressions (|^), (jT^) and (jl^) we can see that all the correlation functions have 
poles at ujq and their residues are expressed as simple combinations of parameters of the 
system. We want to emphasize again that these three results are exact and do not require 
any approximations in their derivations, and they are valid for both the positive and the 
negative U Hubbard model. 

Among all three of these theorems, the third one is the most interesting, since the density 
correlation function can be measured directly in experiments. This theorem predicts that a 
new collective excitation, the r] particle can be found in the density spectrum, whose intensity 
is directly proportional to the square of the superconducting order parameter, in this case, 
an s-wave order parameter. Physically this is because the density correlation function 
measures a particle hole like excitation, and therefore, the t] particle, which has charge two 
rather than zero, can only make a finite contribution if the quasi-particle excitations are 
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a coherent admixture of particles and holes. The mixing amplitude is exactly proportional 
to the square of the superconducting order parameter. For a momentum Q = (tt, tt) the 
energy of the ?7-mode is lvq = U — 2^. However, unlike the first two theorems, the third 
theorem only makes a partial prediction about a correlation function, namely the part that 
is singular at ujq. For this reason, it is highly desirable to find an approximate scheme which 
yields complete information about the full momentum and frequency range of the density 
correlation function, and reduces to the exact theorem for the special momentum Q and 
frequency ujq. Such an investigation is carried out in the subsequent sections. 



III. SELF-CONSISTENT LINEAR RESPONSE THEORY 

A. Formalism 

In calculating the density response in the superconducting state, it is known that one 
must take into account the collective motion of the condensate in a self-consistent manner as 
well as the usual quasi-particle-hole excitations. An externally applied density disturbance 
(f)q at wave length q generates two kinds of responses in the system, single-particle like 
excitations across the energy gap, and a collective motion of the condensate expressed as a 
paring amplitude 77^, and a density field p^, at a finite wave length q. The usual RPA formula 
in the superconducting state takes into account the self-consistent density field p^, h\it docs 
not treat the self-consistent pairing field 7]q. As we shall show, such an approximation is not 
fully self-consistent and violates the SO {A) symmetry of the Hubbard model. On the other 
hand, a fully self-consistent treatment of both the density and the pairing field correctly 
describes the contribution of the 77 particle in the superconducting state, and reduces to the 
exact theorems reviewed in the previous section. 

In the following, we review the basic formalism. In the Appendix we compare our results 
to the formulas one gets from the equations of motion method and from the diagrammatical 
approach. In both cases we find exact agreement. 

Our starting point is the negative-f/ Hubbard Hamiltonian perturbed by an infinitesi- 
mally small external field (j)q which is coupled to the —q component of the density. 

<ij> i i 

We will assume that our system is in a superconducting ground state characterized by the 
nonvanishing averages of the zero momentum operators 

(c.jk^Cfci) = {cl^c^_ki) = UkVk- (14) 
Here the parameters Uk and Vk are given by 
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with efc = —2t{coskx + cos ky) — jj , = y + A'^ and the BCS self-consistency condition 
is imphed 

A = -Uj2up^P (17) 
p 

After we turn on the perturbation, the operators in (14) are no longer the only operators 
with nonvanishing expectation values. The expectation values of the following operators 
carrying momentum q are induced as well. 



ivlit)) 


= E( 

k 




= E( 

k 


(pM) 


= E( 

k 




= E( 

k 



(18) 



The operators in ([T8| ) describe the collective motion of the superconducting condensate. 
Within the self-consistent linear response theory, the local operators respond to the local 
field, which is the sum of the external disturbance (pq and the induced self-consistent field 



given in equation ([TSD. It is very important that in (^) we take into account not only the 
particle-hole channel as is usually done in RPA but also the particle-particle and hole-hole 
channels. 

The Hamiltonian (|T^ can be linearized with respect to the usual BCS pairing fields and 
the induced fields ([T8|) giving : 



i 



per 



+ uip,i + ^) E + u{pq^ + ^) E cI'-.tC'^'t (19) 



^ A:' fc' 



The first line in the last equation will be considered as an unperturbed Hamiltonian Tio- 
Since the fields in (|18|) are proportional to (pq one can use the Kubo formula and treat the 
last four terms in (|T^) as a perturbation Tii. Thus the response (/(t) ) is given by 



(/(t)) = -z/ dt'{ [f{t),n^{t') ) (20) 

and / can be any of the operators ?7^, 77,, or pq^. When combined, these equations form 
a self-consistent set of equations and the expectation values in equation (llSD can be solved 



^ To prevent any possible confusion we want to clarify that in this Section r/g, and rj* = (r/|) 
are the expectation values of the corresponding operators and in commutation relations must be 
considered as c-numbers. 
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for purely in terms of the external field 0g. As an example of how one can proceed from 
(pO|) we will show explicitly the calculations for 77*. 



r]*g{t')dt'Y,{ [ cI+«TWcUiW,C-fc'i(t')cfc'+gT(0 ] ) 

^ kk' 

-'-^ kk' 



°° kk' 



iU 



Mt') + \ct>,{t'))dt'Y.{ [4+,TW4.iW,4'-,T(Oc.'T(0] )^ (21) 

kk' 

When we perform a Fourier transform of the last equation it is helpful to make use of the 
usual trick of replacing the commutator with a time-ordered Green's function and shifting 
the poles to the lower-half plane at the end. This will give us the following expression for 



Vqlo = Jdtriq{t) 



Fourier components of the other functions are defined analogously 



(22) 



where we introduced = Pqcoi+Pqwh Gp{uo) = J e'"^\—i){Tcp^{t)cl^{0))dt and Fp{u) = 
J e^'^^{—i){Tcpi(t)c_pi{0))dt are the usual BCS Greens functions. We remind the reader that 
after performing the frequency integration in ( ^2]) one should shift all the poles of the final 
expression to the lower half-plane. Similar expressions can be easily obtained for rjq^ and 

Pqui- 

Vqu. = -^UY.j ^Fp{u)Fp+q{ + iy + U)7]*q^ 

J2 J ^Gp{v)Gp+q{uj - v)riq^ 

-'UY. J ^Fp{^)Gp+q{y + ^){Pqu> + (23) 

-^UY. [ ^Fp+q{u)Gp{u - 

p J 271 
p J Z7T 

+^UY {Fp{^)Fp+q{^ + ^)- Gp{u)Gp+q{u + UJ)} {pq^ + 0,^) (24) 



Iqu) 



Iqu) 
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In the diagrammatic language these equations correspond to the vertex corrections for 
scattering in the particle-particle, particle-hole and hole-hole channels and are presented in 
Appendix B. Generalization of formulas (p2D - (E4|) to the case of finite temperatures can 



be accomplished by substituting ^ / |^ — ^ Skipping laborious but straightforward 

calculations we present the final expressions leaving their derivation aside. We introduce 
the matrix equation 



V 



-UXbcs 



Vquj ~ VqoJ 
V Pquj + 4>qL0 



(25) 



where the matrix elements are given by the following expressions, 



U+ = UY.[l-f{E, 



■p+q) 



Vtpq{\ 



'"l+q - '"D 



pq 



iO 



+ U[fiE,+,)-fiE,)] 



{Ep+q — Ep){UpVp^q + VpUpj^qY 

- + 20 



p+q) 



fiEp)] 



uj(l 



p+q 



up' - V'^q + ZO 



U\f{Ep^q)~j{Ep)\ 



p+q 



'-O^pq+tO 



UJ 



Uj2[l-f{Ep^q)-f{Ep 



^pqi^^p+q^p ~)~ 'Vp-\-qVp 
Up- - + iO 



P+q) 



fiEp)] 



(^pqi'^p'Vp+q ^p'^p+q) 
- + tO 



f^E [1 - fiEp+q) - f{Ep)] + ^P+'^^P+'l) 



up - v'^^ 



-U\J[E, 



p+q) 



fiEp)] 



^pqiUp-^qVp^q UpVp 



UJ 



pq 



-UY.[l-f{Ep^q)~f{Ep)] 



uiUpVp + Up+qVp+q) 



UJ^ 



u'pq+iO 



+ U[fiE, 



p+q) 



fiEp)] 



^iUp+qVp+q - UpVp) 



uj^-ei 



pq 



iO 



Xbc 



E [1 - fiEp+q) - fiEp)] ^^^^^^^^"^^ + 



Uj"^ - Vpq + «0 



\f{Ep^q) - fiEp)] 



^pqiUpUp-\-q VpVp-\-q) 

uj^ - ei^ + zo 



(26) 



(27) 



(28) 



(29) 



(30) 



(31) 



In the equations above 



pq 



tp+q 



dp , ^pq 



Ep-\-q ~\~ Ep , dpq 



E, 



p+q 



E„ 



and fiE) is the Fermi distribution function. 

Solution of the density part of the matrix equation ( ^5]) can be written in a very suggestive 
form: 
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Pqu) 



-UXbcs + 
f U Xbcs - 



l-t4 



(32) 



-+t++t. 



It is easy to see that the last terms in both the numerator and the denominator are propor- 
tional to the BCS pairing amplitude, therefore, this expression reduces to the usual RPA 
formula in the normal state. In the superconducting state however, these correction terms 
make important contributions, and it would be incorrect to neglect them. These terms take 
into account the self-consistent pairing fields in the superconducting state, and in the di- 
agrammatic language, they arise from multiple scattering in the particle-particle and the 
hole-hole channels. Equation (^2|) agrees with the expression for the density-density corre- 
lation function derived by Kostyrko and Micna^. Formulas (pSj) - ([3^ ) are the core of the 
Self-consistent Linear Response Theory. In Part B we will show analytically that for the 
case oi q = Q and uj = ujq, these formulas reduce to the exact theorem given in the previous 
section, with the correct location of the poles and the residues. In Part C we compare these 
results with numerical calculations. 



B. Comparisons with exact theorems 



Here we show how our formulas (^5]) for T = reduce to the exact results (|10D and (p!2D 
for Q = (tt, tt). This point has a property that QpQ = — 2/i for all momenta p. Therefore QpQ 
may be taken outside the p summation in the equations of (|25|) and a number of interesting 
identities arise. 

We start by multiplying the first equation in ([25| ) by u, the second by flpg and adding 
these two equations. After a few simple manipulations the resulting expression can be 
written as 



vQ 



U{l-n)] ivhui-VQui) 



(33) 



We can prove the following relationship Vpq{vpUp+q + Vp+qUp 
Up+qVp+q) and use it to write the third equation of (^) in the form 

^pqi^p'^p ~l~ '^p+q'^p+q 

P 



2A 



-2A0, 



qijj 



2Ar/ V-^^^-^ 



pq 



Uj"^ - l^pq + iO 



iVquj Vqoj) 



pq 



p ^^-l^^q + iO 



2A 



[UpVp + 



[UpVp + Up+qVp^q) (pqi^ + (j)qi_j) 



(34) 



^ We remind the reader that since the present analysis is for T 
- (^) vanish. 



0, all the fermi functions in (5 
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We add to (p^ the first equation of (pSj) multiplied by QpQ and the second multiplied 
by to obtain^ 

<^ivh. - VQu.) = l-^pQ - U{1 - n)] {7]*Q^ + r]Q^) + (35) 
Equations ( p3D and (|35D can be solved for rig^ and rjg^ in terms of (pg 

where we introduced ujq = U{1 — n) + QpQ = U{1 — n) — 2/i.0 Now, the expressions ( p6D by 
themselves give us a mixed correlation function (|^) or they can be inserted into any of the 
equations of (|25| ) giving the density-density correlation function 

1 + (7XfecHQ,t^) 



at = cuo corresponding to the r^-pair and using 

2A 

we can write (|38|) as 



where hiQ^uj) = J2p "^+^^"^+'3 . From ( p8D we immediately see that D{Q,u) has a pole 



1 + Uxbcs{Q,u;) = - — - ^po)^2(g,^) (39) 



A2 1 

D(Q,uj) = 4:-—r- \- part regular at uo (40) 

U^{l — n) u — Uq + zO 

Thus the Self-consistent Linear Response expressions agree with the exact theorems, giving 
the correct position of the resonance as well as the same values for the residues of all 
correlation functions. Here we have explicitly verified theorems 2 and 3 from last section. 
In order to verify theorem 1, one needs to add an external pairing disturbance, in addition 



^ Equations (|3|) and 

and J2p{Vpq + Vpq) 

^ This expression may seem to give a different value for than @ . This comes from the fact that 
in the Linear Response we do not take into account the Hatree-Fock corrections to the self-energy 
or in other words in the diagrams of Appendix C we do not consider corrections to the single 
particle Green's functions. In fact, the only effect of these corrections is to renormalize a chemical 
potential by an average effective field on each particle due to the interaction with the particles 
of opposite spin fi — > ij, — After such substitution this last expression for ujq is exactly the 
same as in (^. This procedure of Hartree-Fock corrections reducing to the renormalization of /i is 
explicit in the equations of motion. 



) arise directly if one writes the equations of motion for J^piVpq ~ Vpq 
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to the density disturbance. Repeating the manipulations given above, it can be seen that 
theorem 1 is verified as welL 

These resuhs show that while the Self-consistent Linear Response Theory is still an 
approximate scheme, it is conserving in the sense that it respects the S'0(4) Wards identities 
expressed through theorems 1 to 3. On the other hand, a simple RPA expression in which 
one neglects the last terms in the numerator and the denominator of equation (|3^) violates 
the 5*0 (4) symmetry and is a poor approximate scheme. 



C. Numerical results 

Here we give numerical results for the spectral functions calculated from equation (p2|), 
with the full momentum, energy and temperature dependence. The calculations have been 
carried out for different t/s' and n's and the values of parameters are specified in the captions 
for each figure. The imaginary infinitesimal in the denominators was taken to be F = 
1.0 X 10~^ and the integrations were performed on a 1024 x 1024 lattice. 

In figure (|I|) we plot the imaginary part of the density response function, namely the 
spectral function, at zero temperature and q = Q. In curve (a) we show the spectral function 
calculated from the Self-consistent Linear Response theory, and show that the location of 
the peak agrees exactly with the value predicted by the exact theorem, indicated by a delta 
function arrow. The width of the peak is not intrinsic, but due to the finiteness of F used 
to smooth the calculations. In curve (b) we show the results from the usual RPA formula, 
and we see that it completely misses the location of the collective mode. For a larger value 
of U, this discrepancy becomes even stronger. In curve (c) we show the result of the simple 
BCS quasiparticle approximation to the density response function. It has an onset at the 
minimal value of Ep + -Ep+g = 2|/i|, but no sharp peak structure of the collective mode. 

In figures (||) and @ we show the momentum dependence of the peak, by plotting the 
spectral function for various values of the center of mass momentum. It is interesting to 
note the difference in behaviour of the peak for different values of the interaction strength 
U. For the case of smaller ?7's shown in Figure we see that the collective mode is most 
sharply defined at Q = (tt, vr), and that it broadens and disperses downwards in energy as 
one moves away from Q. For larger f/'s, however, the energy of the mode increases as we 
go away from (vr, vr) with only a small change in the width or intensity of the peak. Within 
a simple T-matrix approximation of the 77-mode, i.e. without the mixing into the density 
excitations, one can show that the ?7-mode always disperses downwards as one moves away 
from the (7r,7r) point, and it merges tangentially into the particle-particle continuum. The 
upward dispersion of the ?7-mode for large U can be attributed to the strong mixing of the 
?7-mode with density excitations. We present a more detailed discussion of the dispersion of 
the mode and comparison between T-matrix approxiamtion and the Self-consistent Linear 
Response method in Section |V|. 

It is important to remind the reader that the special properties of the momentum Q 
demonstarted in figures above have nothing to do with the nesting property of the fermi 
surface, since we are studying a doped system where 2k f is markably different from Q at 
this filling. Instead, the special role of Q is that at this momentum, the rig operator becomes 
an eigenoperator of the Hamiltonian. Physically speaking, this is the momentum where the 
particle-particle continuum vanishes. 
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In figure (^) we plot the temperature dependence of tlie spectral function for momentum 
Q. We see that the peak intensity decreases with temperature and the peak vanishes exactly 
at Tf.. This behavior is in exact agreement with the third theorem that the peak intensity is 
proportional to the BCS order parameter. We also see that at low temperatures, the weight 
of the collective mode is transferred from the higher energy continuum. 
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FIG. 1. ImD{Q,uj) versus oj calculated using BCS, Random Phase Approximation and 
Self-consistent Linear Response Theory formulas. The arrow indicates ojq - s, position of the jy-pair 
as predicted by the exact theorems. This plot is for U = —t, n = 0.79 and T = 0. 
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FIG. 2. ImD{q, u) versus a; for ?7 = —t, n = 0.87, T = and different momenta q. The arrow 
marks lvq. This figure shows that for small C/'s the collective 77- mode exists only in a very close 
vicinity of (tt, vr). One can also see a small dispersion of the mode - as we go away from (tt, tt) the 
energy decreases. 
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FIG. 3. ImD{q,oj) versus u for U = —4t, n = 0.87, T = and different momenta q. The 
arrow marks ujq. Notice a big positive dispersion and absense of strong damping of the mode away 
from (tt, tt). 
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FIG. 4. ImD{Q,uj) versus lu for U = —t, n = 0.87 and different temperatures: T = 
corresponds to A = 0.041t; T = O.OlSt to A = 0.031t; T = 0M7t to A = 0.021t; T = 0.019t to 
A = 0.012t and T = 0.0197t to A = 0.002t. The arrow marks loq. It is interesting that as we raise 
the temperature, the intensity of the peak decreases but its position does not change. 
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IV. MONTE CARLO CALCULATIONS 



It is useful to compare the results of the exact theorems and the self-consistent cal- 
culations with the numerical data from Quantum Monte Carlo simulations. Using the 
finite-temperature determinantal Monte Carlo technique0 we have measured the Matsubara 
time-ordered Green's functions 

D{q,T)=<Trp.,{T)p,{0)> (41) 

and 

Piq,T)=<T,r^^{T)r^l{0)> . (42) 

The spectral weights lmD{q,uj) and lmP{q,uj) can be obtained by maximum entropy ana- 
lytic continuationEJ'Ej of D{q, r) and P{q, r) to the real frequency axis. 

In this section, we will present results for U = — 4t and —It on an 8 x 8 lattice. We 
will consider electron fillings (n) = 0.87, 0.7 and 0.5. For a filling of 0.87 and U = — 4t, 
the Kosterlitz-Thouless superconducting transition temperature T^^ is ~ In this 

parameter regime, we have carried out calculations for temperatures varying from T = 0.5i: 
down to 0.08t. 

The solid line in Fig. (^) shows Im P{q, u) at center-of-mass momentum q = (vr, vr) for 
U = —it, (n) = 0.87 and T = 0.33t. Here the chemical potential is /i = — 2.15t. We see 
that the peak is centered at ~ 0.37t, which is close to U — 2/i. The finite width of the 
peak is due to the resolution of the maximum entropy method. The dotted and the dashed 
curves in Fig. (|^) show lmP{q,uj) for q = (7r,37r/4) and (37r/4, 37r/4). We observe that 
the 7] resonance broadens and shifts to higher frequencies as q moves away from (tt, vr). The 
excitation is expected to disperse as \q — (vr, vr)p, however, the momentum resolution on the 
8x8 lattice is 7r/4, and hence these calculations can't provide information on the dispersion 
of the ?7-resonance in the immediate vicinity of (vr, vr). 

Fig. (H) shows the spectral weight of the density response function, lmD{q,Lj), at 
q = (vr, vr) for various values of T. As the temperature is lowered a sharp peak developes 
at low frequencies. We note that the position of the peak in lmD{q = {ii,7i),uj) coincides 
with the position of the peak in ImP(g = {it , it) , u) . In the grand canonical ensemble, the 
chemical potential depends on the temperature, and this is reflected in a slight shift of the 
peak position. 

We can obtain further information on the density fluctations by studying D{q,uj = 0). 
Fig. (^ shows D{q,0) versus q for q around the Brillouin zone at various temperatures. 
According to the Kramers-Kronig relation 

,~<iyimD(^_ (43) 

JO vr uj' 

D{q, 0) provides us with information on the integral of the spectral weight divided by u. In 
this figure we see that D{q = (vr, 7r),0) increases rapidly as T is lowered. We believe that 
this rapid growth of D{q = (vr, vr), 0) is due to the coupling of the density excitations to the 
T] particle-particle channel as T — T^^. However, for q ^ (vr, vr) such a rapid growth is 
not observed. This would then mean that the ri excitation is relatively local in momentum 
space. 
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In Figure (j^), we have seen that the peak in ImD starts growing for temperatures 
considerably higher than T^"^. This is possibly due to the fact that there are significant 
superconducting fluctuations above T^f". In order to obtain information on the strength of 
these superconducting fluctuations, we consider the current-current correlation function 

1 ri3 

Axx(g, iu^m) = ^ / rfre*"™" < Jx(g, r)j,{-q, 0) > (44) 
iV JO 

with 

Jx {q) =itJ2 e'^^'^icl+xaCea - cLq+x.) • (45) 

e 

The superfluid density Ds is given byil 

^ = -{kx) - ^xx{qx = 0,qy ^ 0, iuJrn = O) (46) 
ne 

where (kx) is the average kinetic energy in the x direction. In Fig. (||) we have plotted 
— {kx) —-^xx{qx = 0, Qy, ioJm = 0) versus qnat various temperatures. These results, which are 
identical to those given in Fig. (8) of Ref.E3, show that there are significant superconducting 
fluctuations on the 8x8 system for T higher than T^"^ ^ O.lt. Basically the coherence 
length of the superconducting fluctuations is extending over the finite 8x8 lattice. We also 
note that the growth of Dg and the peak in ImD appear to be strongly correlated, as one 
would have expected from the exact theorems of Section II. 

For U = —4t the statistical error in the Quantum Monte Carlo data is large, and as a 
result of this the spectral weights get smeared by the maximum entropy analytic continu- 
ation. For instance, in Fig. we have seen that the t] resonance, which is a 5 function, 
was broadened considerably, and its frequency was shifted. However, for U = —It we can 
obtain Monte Carlo data with better statistics and the resulting spectral weights have im- 
proved resolution. Unfortunately, for U = —It, the superconducting transition takes place 
at temperatures much lower than it is feasible for a numerical study. Hence, we will only 
study the r] pair field susceptibility at the elevated temperature of T = 0.25t. 

Here we consider the filling dependence of the r] resonance. Fig. (|^) shows lmP{q = 
{n,n),uj) versus u/t for fillings of 0.87, 0.7 and 0.5. The corresponding chemical potentials 
are — 2.16t, — 2.67t and — 3.25t. We see from this figure that the resonance frequencies are 
very close to the expected values of U — 2/i. In addition, we find that the total spectral 
weights in the peaks are 0.14, 0.31 and 0.49, respectively. These values agree well with the 
exact theorem of Eq. (8), which gives 1 — (n) for the weight. 

In this section, we have seen that a sharp resonance exists in the r^-pairing channel. As 
the temperature is lowered and the superconducting correlations develop across the 8x8 
lattice, the density fluctuations can couple to this channel. The resulting density excitation 
spectrum has a peak at the resonance frequency in addition to a continuum of excitations. 
These results are in agreement with the predictions of the exact theorems discussed in Section 
II and the results of the self-consistent calculations presented in Section HI. 
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FIG. 9. ImP((7 = (vrjVr),^^) versus ujjt for U = —It, T = 0.25t and various values of (n) 
V. MAPPING OF THE ??-MODE TO ANTIFERROMAGNETIC SPIN WAVE 



There is an elegant way to derive the dispersion of the 77-mode for the negative U Hubbard 
model in the large \U\ limit by mapping it onto a Heisenberg model in an external mag- 
netic field. Mele first used this formalism to discuss the dispersion of the r] mode for zero 
temperatureii. Below we generalize this approach to include the effects of finite temperature 
and obtain the temperature dependence of the rj mode dispersion. 

Our starting point is a particle-hole transformationBUi^ on the bipartite lattice: 



where A and B stands for the two sublattices. In the future we will write the transformation 
of the spin-down particles in the form Cj| (— )*c||. 
For the two-fermion operators (^7]) gives 



A. Particle-hole transformation 




(47) 



(48) 



and so the charge Q = + rii becomes the spin Sz = |(^t — ni) and vice versa. 
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Q^2S, + 1 2S,^Q-1 (49) 

If we now write our hamiltonian in the form 

^ = E 4aCja - \U\ E(^iT - hina -h-f^Yl ^icia (50) 

<ij>cT i ia 

which is different from the earher notation by only a chemical potential renormalization 
jl = fx + we obtain using (^81) 

<ij>cr i ia 

Hamiltonian (^) describes a positive U Hubbard model at half-filling [| in the presence of 
a magnetic field 2fi Q which in the large U limit is known to be equivalent to the Heisenberg 
model 

n = JYl s,-Sj -hY^Si, (52) 

<ij> i 

with J = ^ and h = 2jl. Therefore by considering the Heisenberg model, we will be able 
to get a deeper insight into the modes of the negative U Hubbard model. 

Our next step will be to establish the correspondence between excitations of the two 
systems. To do so, we notice that the particle-hole transformation turns SU{2)pseudospin 
algebra into a usual SU{2)spin algebra 



= 


i 




= E^ltCii 

i 


V = 


E(-)'caCiT 

i 


s. 


= E 4iCiT 

i 


Vo = 


-\{Ne-N) 




i 



(53) 

And it transforms Cooper pairs and charge density wave into antiferromagnetic waves in 
X — y and z directions respectively. 

= E 4Ai ^ s^iQ) = T.H%c,i 

i i 
i i 

P{Q) = E (4tC.T + - 2S,iQ) = Yi'Y (c1tC.T - 4lC^l) (54) 



^ We refer the reader to a well known theorem that /u = ^ at half-filling 
* Here and everywhere else we take Bohr's magneton to be one, /i^ = 1. 
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Hamiltonian (|52D, as well as any other rotationally invariant system subject to an external 
magnetic field, has a Larmor precession mode given by 

[n, S±] = ±hS± (55) 

So, from ( |53[ ) we see that 77-excitation is what becomes of the Larmor motion under the 
particle-hole transformation. 

Coupling of the 77 excitation to the density fiuctuations also has its analogy in the Lar- 
mor motion. We assumed s-wave superconductivity in the negative U case, then after the 
particle-hole transformation takes us to the Heisenberg model we must have an antiferromag- 
netic order in x — y plane or non- vanishing < S±{Q) >. The BCS order parameter being real 
implies that antiferromagnetism will be in x-direction. This means that [S±, Sz{Q)] = S±{Q) 
is a c-number and so the Larmor motion is conjugate to the Sz{Q) oscillations. Thus when 
we excite Larmor oscillations we will have a response in S^iQ)- This is analogous to the 
negative U model in which rj and CDW modes are conjugate 7]\ p{Q) = and therefore 
coupled to each other. 

From the fact that the {^7^ ''7, p((5)} operators for the negative U Hubbard model corre- 
spond to {S+, S^, Sz{Q)} for the Heisenberg model, it follows that oscillatons with wavevec- 
tors close to (vr, vr) will also be in one-to-one correspondence for the two systems and will 
therefore have the same dispersion. So in the next section we briefiy review what is known re- 
garding the excitations of the spin-wave modes of the Heisenberg Hamiltonian in an external 
magnetic field, Eq. (^). 

B. Dispersion for the Heisenberg model 

We consider a mean field spin-fiop groundstate of the Hamiltonian, Eq. (|5^ , which 
has antiferromagnetic order in the x-direction and is polarized in the 2;-direction by the 
magnetic field h so that 

<S.> =-l 1. (56) 

As previously discussed, the antiferromagnetic order corresponds to the superconducting or- 
der in the negative U system and the magnetic polarization to the deviation of the band filling 
from half-filling. This means that temperature-wise <y{T) oc A(r) and /3 oc 1 — n = const. 
If we now write the equations of motion for the individual spins 

^ = m,S,] (57) 
and linearize them around < Sj > we get 

at 2 

^ <j> ^ <j> 



dt 



dt ^ ' 2 
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where < j > stands for the sites nearest to i. We now take Mf 
and perform a Fourier transformation of ( ^8]) to obtain 



'5S' and M 



J[3 



(4 - lk)Ml + hMt 



In 

-^^ML = -^(4 + 7.)ML 



(59) 



with = X^a c**^- Equation (|59D may be solved giving the dispersion relation^ 



2^,2 



;i6-7fc) + 



(4 - 7fc) - 



For small /c's the last expression may be expanded to give 

ujl = h'^ + {2J^a^ - hp J) \{Kay + {kyo) 



(60) 



(61) 



Now we are in a position to discuss dispersion of the r^-mode and compare the Self- 
consistent Linear Response with the T-matrix approximations. Above Tc, there is no mixing 
between the particle hole and the particle particle channels, a simple T-matrix calculation 
yields the correct answer of a downward dispersion. This result agrees with equation (BTl) 



when we take oc |Ap = above T^.. In the superconducting state a simple T-matrix 
approximation without taking into account the mixing is no longer valid. In fact, we see 
from ( |6lD that the mixing term U'^o? gives a positive contribution to the dispersion. These 
two different terms compete with each other, so that a positive dispersion can be obtained 
at zero temperature for large |?7|. On the other hand, for small |[/|, the second term 
—hpj always wins over the first term 2J'^o? and the dispersion is negative even at zero 
temperature. This is consistent with the numerical results of the Self-consistent Linear 
Response theory, presented in Section [111 L] . From (pT]), we also see that at large |f/|, as one 
rises the temperature from zero to Tc the temperature dependence of a leads to a continious 
transition from the positive dispersion to negative. 



VI. SUMMARY AND CONCLUSIONS 

In this section we would like to summarize the results obtained in the previous sections, 
and comment on their implications. We see that the analytical approaches based on the 
50 (4) symmetry property of the Hubbard model, the self-consistent linear response the- 
ory and the numerical Monte Carlo simulations agree with each other. These calculations 
establish that the ri particle has the following properties: 

1) It has charge quantum number two, and is a sharp excitation of the system for a 
finite range of total momentum around Q = (7r,7r), and the mode disperses as the total 
momentum departs from Q. The special role played by Q is independent of the nesting 
property of the fermi surface. 

2) It has spin zero and energy U — 2fi. 
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3) Because it has charge two, it contributes to the particle-hole fluctuation spectrum 
(in this case the density spectrum) only when the superconducting pairing correlations are 
present. The analytic calculations show that the intensity of the peak onsets as the square 
of the superconducting order parameter. The energy of the peak is unchanged as a function 
of temperature. 

The value of present detailed study of the rj particle is to use it to illustrate a mechanism 
with which particle particle excitations can couple to the particle hole spectrum below Tc, 
and to test the methodology of the theoretical approach. We see that the Self- consistent 
Linear Response method correctly takes into account the contribution of the particle-particle 
excitation to the density fluctuation spectrum, and agrees with both the exact theorems 
derived from the 5*0 (4) Wards identities of the Hubbard model and the numerical Monte 
Carlo simulation results. This method is general, and thus provides a useful starting point 
for investigating the collective modes in the superconducting state. 

In the next paper, we shall use the methodology developed and tested in this paper to 
carry out an investigation of the tt particle of the positve U Hubbard or the t — J model, 
which has spin one rather than zero, but shares many other properties of the r) particle. 
We shall discuss in detail the similarities and differences of these two collective modes and 
compare our results with the neutron scattering experiments. 
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VIII. APPENDIX A: COMPARISON WITH EQUATIONS OF MOTION 

METHOD 



In the equations of motion method we introduce 



Vkq — cl+q|clfc| 
Vkq = C-k-qlCk'\ 

Pkq = Pkql + Pkql = Cfc+gt'^'^T + ^^-kl^^-k-ql 



(62) 



and write the Heisenberg equations of motion for these operators using (|T^) as the Hamilto- 
nian. These equations of motion are then factorized in terms of the occupation numbers for 
the electrons and BCS anomalous averages. The details of such computations are presented 
extensively in the literaturefr^'i and so we will only outline the main steps without going 
into details or discussing physical meaning of the equations. 



[n,Ppg^] = UJpqPpq^ - {vl^g - vl)U ^pkqi + ^^gj 

+ UpVpU Viq - Up+qVp+gU Vkq + ^{vl 
k k 

[H, Ppql] = -^pqPpqi + (^^p+g " vl)U | ^ PfcgT + ^</'<?| 



~%q) 



+ Up^qVp+gU viq - UpVpU Y Vkq + HV^q " Vpq) 



(63) 



(64) 



n,vU =npqvU + u{i 



2 1 
^p+q - 



J2Vkq 



+ UUp+qVp+g l^PfcgT + \<Pq^ + UUpVp | ^ Pfcgi + ^(pq^^ + ^{ppq^ + Ppql) (65) 



[n, ripq] = -npqr]pq - U {1 ~ v^^^ - f J) ^ r]kq 



- UUp+qVp+q jx^Pfcgt + \(t^}^ - UUpVp jx^Pfcgi + ^'/'gj " ^(PpgT + Ppqi) (66) 

And we have already included the Hartree-Fock terms into the renormalized chemical 
potential. Following Anderson we take the second-order commutators 



'H,[n,npq-nlq]\ ={^l^ + A/\'){n_ 



pq Vpq 

+ UVLr,q{l - vl^^ - vl) YiVkq - vi 



pq 
pq\ 



2AuJpq{ppq1 - Ppql) 



+ 2UA{UpVp + Up+qVp+q) YiVkq " Vkq 



- 2UA{n 



p+q 



Y(Pkq1 ~ Pkql) 



U [UpVp + Up+qVp+q) 



'^i^iPkql + Pkql) 
k 
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-U{l-u 



2 _ 2^ 

p+q ~ "^p) 



Y^iVkq + via 



T-C, [H, Vpq + Vlq]] = ^IqiVpq + Vlq) + 2AQpg{ppg^ + Ppg^) 

+ Uflpg{l - Vl^g - J2iVkq + Viq) 

k 

+ UVLpq{UpVp + Up+qVp+q) < Y^{pkq\ + Pkqi) + 



- U{1 - Vlq - U 



Y.(^kq - Vlq 



[n, [H,ppq^ - ppq^ = UJl 



- Ppql) - 2AuJpq{ripq - T^Jj 

- UuJpqiUpVp + Up+qVp+q) ^{r]kq " 

k 

+ UUpqiyl^^ - Y^iPkq] - pkqi) 



U{v'q - vl) 



'^i^i.Pkq] + Pkqi) 



+ U {UpVp — Up^qVp^q) 



^^Y.^Vkq + vlq 



[H, [n, Ppq^ + Ppql]] = + 4:A'^){ppqi + Ppqi) + 2Afipg(r7pg + T], 



pqJ 



+ UuJpq{UpVp - Up+qVp+q) Y^{r]kq + 



+ 2AU{1 - V. 



2 

p+q p 



J2(Vkq + Viq 



+ 2AU{UpVp + Up+qVp+q) |Xl(PfcgT + Pkqi) + 0g | 



+ U{v. 



p+q "^p) 



'^i^i.Pkqi - Pkqi) 



U {UpVp + Upj^qVpj^q) 



KY.(Vkq-viq 



(67) 



(68) 



(69) 



(70) 



Now we take [H, f] = —ujf, where / can be any of p, t]'^ or b, combine equations (|67D 

Op ) cancel out and get after 



with (|69|) and (^) with (|70D such that unphysical poles (cj^ 
a few straightforward transformations: 

(^^ - ^lq)Ppg = UVLpqiUpVp + Up+qVp+q) Y^iVkq + Vkq) 

k 

- Uuj{UpVp + Up+qVp+q) Y^iVkq " Vkq) 

k 

+ U Vpq{VpUp+q + Vp+qUp)'^(Y^ Pkq + 0g) 



(71) 



(^^ - ^Iq) ivlq + Vpq) = UQpqil 



vl+q-vl)Y,{viq+Vkq) 
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-UUJ{1- Vl^g - VD Y^ivlg - Vkq) 

k 



+ U ^Ipq {UpVp + Up+qVp+q) Pkq + (f)q) (72) 

k 

k 

+ U Upq{UpUp+q + VpVp+qf J^iVkq " Vkq) 

k 

- Uuj{UpVp + Up+qVp+q)(^ Pkq + 0g) (73) 

k 

One can see that after we divide equations ( [71]) - (^) by (u;^ — i/^^) and sum over p we recover 
the formulas of the Linear Response. The equations of motion method is very often more 
convenient because one gets equations before they were summed over momentum and this 
gives some additional freedom in manipulating them. So, for example it is much easier to 
see the //-pair energy-eigenvalue or to compute the residue of -D(g, oj) at ujq in equations(|7T]) 



- ([73|) than in (pSj) - (0). The same time the Linear Response Theory has an advantage of 



being more intuitively clear and transparent. 
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IX. APPENDIX B: DIAGRAMMATIC EQUATIONS 



In this appendix we show the diagrams that correspond to equations (12^) - 
The first equation describes scattering in the particle-particle channel and should be 



compared to (|22|) 





+ 




+ 




+ 




This next equation is for scattering in the hole-hole channel and corresponds to 





+ 




+ 




+ 
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The two following equations are for scattering in the particle-hole channels with 
or down. Their sum gives (24). 



Spin up: 




T F ^ 



Til T 



— . T 



Til T 



i G ^ i G ^ T F t 



Spin down: 




i\ i F ^ 



F* G F* 

^ G G If 
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